function City_Plots_20231030

  % Assume index is observed at the end of the file month.
  %
  % This program:
  %
  % 1. Computes summary statistics for the city price, rent and price/rent data
  %
  % 2. Plots Real Price and Real Rent for individual cities
  % 
  % 3. Plots Real Price and Real Rent for groups of individual cities
  %
  
  tic
  
  computer = 2;
  if computer == 0 % mac
    addpath(genpath('/users/fefama/Dropbox/Z/Papers/Real Estate'))
  elseif computer == 1 % pc
    addpath(genpath('C:\Users\fefama\Dropbox\Z\Papers\Real Estate'))
  elseif computer == 2 % pc
    addpath('d:\dropbox\projects\Real Estate', 'd:\dropbox\projects\Real Estate\matlab', 'd:\dropbox\projects\Real Estate\Data 20230103' )
  end

  c.Oct_2022 = 430;                % 202310
  p.last_mth = c.Oct_2022;         % Sep_2019 or Oct_2020 or May_2021

  c = setup_global_constants(c, p);
  p = setup_parameters(c, p);
  
  p.SA_prices = true;
  p.SA_rents  = true;     % Relevant only for All

  p.CPI_SA    = true;     % true: real prices and rents  divided by CPI_SA


  % p.fid is output file
  % 
  [price_filename, rent_filename] = get_input_file_names(c, p);


  % load data
  %
  d = load_CPI_price_and_rent_data(price_filename, rent_filename, c, p);
  

  % setup output file...p.fid is output file
  %
  p = setup_output_file_etc(c, d, p);
   

  % write raw data
  %
  if p.wrt_raw_data
    wrt_raw_data(d, c, p);
  end

  % Plot Figure 1
  %
  city = 1;
  cntrl_plot_city_real(city, c, d, p);
  toc
end

%-----------------------------------------------------------------------------------------
%                                                                                        *
%                       Functions to set up constants, load data, etc.                   *
%                                                                                        *
% ****************************************************************************************
%
function c = setup_global_constants(c_in, p)

  c = c_in;
  
  c.num_2   = 2;
  c.num_3   = 3;
  c.num_4   = 4;
  c.num_5   = 5;
  c.missing = -99.99;

  c.all     = 1;
  
  c.price   = 1;
  c.rent    = 2;
  c.ratio   = 3;

  c.nom     = 1;
  c.real    = 2;
  
  c.price_to_rent = 1;
  c.rent_to_price = 2;


  % Assume index is observed at the end of the file month.
  %
  c.Jan_1987 = 1;     % 198701
  c.Jan_1990 = 37;
  c.Jan_1991 = 49;    % 199101
  c.Dec_1996 = 120;
  c.Dec_1997 = 132;
  c.Jan_2000 = 157;   % 200001
  c.Feb_2006 = 230;   % 200604
  c.Mar_2006 = 231;   % 200604
  c.Apr_2006 = 232;   % 200604
  c.May_2006 = 233;   % 200605
  c.Mar_2009 = 267;   % 200904
  c.Apr_2009 = 268;   % 200904
  c.May_2009 = 269;   % 200905
  c.Jan_2012 = 301;   % 201204
  c.Feb_2012 = 302;   % 201204
  c.Mar_2012 = 303;   % 201204
  c.Apr_2012 = 304;   % 201204
  c.May_2012 = 305;   % 201205
  c.May_2019 = 389;   % 201905
  c.Sep_2019 = 393;   % 201909
  c.Apr_2020 = 400;   % 202004
  c.Oct_2020 = 406;   % 202010
  c.Nov_2020 = 407;   % 202011
  c.Dec_2020 = 408;   % 202012

  c.frst_mth     = c.Jan_1987;
  c.last_mth     = p.last_mth;
  c.num_mths_ttl = c.last_mth;
  c.base_date    = 198701;
  
    c.real_key_obs = [c.Jan_1987; c.Dec_1997; c.Jan_2000; c.Mar_2006; c.Apr_2009; c.Feb_2012; c.last_mth];

  c.rescale_1987_last_mth_ave =  0;
  c.rescale_frst_5_yr_ave     = -1;
  c.rescale_frst_10_yr_ave    = -2;

  c.USAve = 1;  c.BMA = 2;  c.NYNY = 3;  c.MFL = 4;   c.DMI = 5;  c.CIL = 6;  c.AGA = 7;  c.DTX = 8;  c.SWA = 9;  c.SFCA = 10;  c.LACA = 11;   % cities with both prices and rents
  c.SDCA = 12;  c.PAZ = 13; c.DCO = 14; c.LVNV = 15; c.WDC = 16; c.POR = 17; c.CNC = 18; c.COH = 19;  c.TFL = 20;  c.MMN = 21;                 % cities with just prices

  
  

  % c.cityabrev and codes for cities (e.g., c.USAve = 1) must be linked one-for-one
  %
  % c.order_of_cities_stored is the order the data are stored; Cities with both prices and rents MUST be first; Cities with rents but not prices are not loaded.
  %                               1        2      3       4      5      6      7      8      9      10      11      12      13     14     15      16     17     18      19      20      21
  c.order_of_cities_stored   = [c.USAve, c.BMA, c.NYNY, c.MFL, c.DMI, c.CIL, c.AGA, c.DTX, c.SWA, c.SFCA, c.LACA, c.SDCA, c.PAZ, c.DCO, c.LVNV, c.WDC, c.POR, c.CNC,  c.COH,  c.TFL,  c.MMN];  % include only cities that are used
  c.cityabrev                = ['USAve'; 'B_MA '; 'NY_NY'; 'M_FL '; 'D_MI '; 'C_IL '; 'A_GA '; 'D_TX '; 'S_WA '; 'SF_CA'; 'LA_CA'; 'SD_CA'; 'P_AZ '; 'D_CO '; 'LV_NV'; 'W_DC '; 'P_OR '; 'C_NC '; 'C_OH '; 'T_FL '; 'M_MN '];


  c.order_for_output_price = [c.USAve, c.BMA, c.NYNY, c.WDC, c.TFL, c.MFL, c.MMN, c.DMI, c.CIL, c.COH,  c.DCO,  c.LVNV, c.CNC, c.AGA,  c.PAZ, c.DTX, c.SWA,  c.POR, c.SFCA, c.LACA, c.SDCA]; % include only cities that are used
  
  c.num_cities_price   = 21;
  c.num_cities_both    = 11;


  % order_of_cities_stored tells the order of cities in the data...which city comes next...from location to city...Cities with both prices and rents MUST be first
  % order_for_output tells the order of cities in the output...which city comes next...from location to city
  % city_loc_in_data tells where a city is in the data...which column...from city to location
  % city_loc_in_output tells where a city is in the output...which row...from city to location

  c.city_loc_in_data         = zeros(1, c.num_cities_price);
  c.city_loc_in_output_price = zeros(1, c.num_cities_price);
  for city = 1:c.num_cities_price
    c.city_loc_in_data        (c.order_of_cities_stored(city)) = city;  % where is a city in the data...which column for a city
  end
  
  for city = 1:c.num_cities_price
    c.city_loc_in_output_price(c.order_for_output_price(city)) = city;  % where is a city in the output...which row
  end
  

  %               1        2        3        4        5        6         7        8        9       10       11       12       13       14       15       16       17       18       19       20       21
  c.name_5  = ['USAve'; 'B_MA '; 'NY_NY'; 'M_FL '; 'D_MI '; 'C_IL ';  'A_GA '; 'D_TX '; 'S_WA '; 'SF_CA'; 'LA_CA'; 'SD_CA'; 'P_AZ '; 'D_CO '; 'LV_NV'; 'W_DC '; 'P_OR '; 'C_NC '; 'C_OH '; 'T_FL '; 'M_MN '];

  c.name_7  = ['USA Ave'; 'Boston '; 'NY     '; 'Miami  '; 'Detroit'; 'Chicago'; 'Atlanta'; 'Dallas '; 'Seattle'; 'SF     '; 'LA     '; ...
               'SD     '; 'Phoenix'; 'Denver '; 'LV     '; 'DC     '; 'Portlnd'; 'Chrlott'; 'Clvland'; 'Tampa  '; 'Minneap'];

  c.name_13 = ['USA Ave      '; 'Boston       '; 'New York     '; '        Miami'; 'Detroit      '; 'Chicago      '; 'Atlanta      '; 'Dallas       '; 'Seattle      '; 'San Francisco'; 'Los Angeles  '; ...
               'San Diego    '; 'Phoenix      '; 'Denver       '; 'Las Vegas    '; 'Washington DC'; 'Portland     '; 'Charlotte    '; 'Cleveland    '; 'Tampa        '; 'Minneapolis  '];



  % c.full_city_name must line up with c.name_5 and codes for cities (e.g., c.USAve = 1) must be linked one-for-one
  %
  c.full_city_name = strings(c.num_cities_price, 1);
  
  c.full_city_name(c.USAve, 1) = "US Average";
  c.full_city_name(c.BMA,   1) = "Boston";
  c.full_city_name(c.NYNY,  1) = "New York";
  c.full_city_name(c.MFL,   1) = "Miami";
  c.full_city_name(c.DMI,   1) = "Detroit";
  c.full_city_name(c.CIL,   1) = "Chicago";
  c.full_city_name(c.AGA,   1) = "Atlanta";
  c.full_city_name(c.DTX,   1) = "Dallas";
  c.full_city_name(c.SWA,   1) = "Seattle";
  c.full_city_name(c.SFCA,  1) = "San Francisco";
  c.full_city_name(c.LACA,  1) = "Los Angeles";
  c.full_city_name(c.SDCA,  1) = "San Diego";
  c.full_city_name(c.PAZ,   1) = "Phoenix";
  c.full_city_name(c.DCO,   1) = "Denver";
  c.full_city_name(c.LVNV,  1) = "Las Vegas";
  c.full_city_name(c.WDC,   1) = "Washington";
  c.full_city_name(c.POR,   1) = "Portland";
  c.full_city_name(c.CNC,   1) = "Charlotte";
  c.full_city_name(c.COH,   1) = "Cleveland";
  c.full_city_name(c.TFL,   1) = "Tampa";
  c.full_city_name(c.MMN,   1) = "Minneapolis";

  c.full_city_name_char(c.USAve, :) = 'US Average   ';
  c.full_city_name_char(c.BMA,   :) = 'Boston       ';
  c.full_city_name_char(c.NYNY,  :) = 'New York     ';
  c.full_city_name_char(c.MFL,   :) = 'Miami        ';
  c.full_city_name_char(c.DMI,   :) = 'Detroit      ';
  c.full_city_name_char(c.CIL,   :) = 'Chicago      ';
  c.full_city_name_char(c.AGA,   :) = 'Atlanta      ';
  c.full_city_name_char(c.DTX,   :) = 'Dallas       ';
  c.full_city_name_char(c.SWA,   :) = 'Seattle      ';
  c.full_city_name_char(c.SFCA,  :) = 'San Francisco';
  c.full_city_name_char(c.LACA,  :) = 'Los Angeles  ';
  c.full_city_name_char(c.SDCA,  :) = 'San Diego    ';
  c.full_city_name_char(c.PAZ,   :) = 'Phoenix      ';
  c.full_city_name_char(c.DCO,   :) = 'Denver       ';
  c.full_city_name_char(c.LVNV,  :) = 'Las Vegas    ';
  c.full_city_name_char(c.WDC,   :) = 'Washington   ';
  c.full_city_name_char(c.POR,   :) = 'Portland     ';
  c.full_city_name_char(c.CNC,   :) = 'Charlotte    ';
  c.full_city_name_char(c.COH,   :) = 'Cleveland    ';
  c.full_city_name_char(c.TFL,   :) = 'Tampa        ';
  c.full_city_name_char(c.MMN,   :) = 'Minneapolis  ';

  c.full_city_name_label = c.full_city_name;
  c.full_city_name_label(1, 1) = "{\itAll}";
  

  % build header with city abbreviations
  %
  c.header_cities = reshape(blanks(c.num_cities_price * 8), c.num_cities_price, 8);
  for row = 1:c.num_cities_price
    c.header_cities(row, :) = sprintf('   %5s', strjust(c.name_5(row, :), 'right'));
  end      


end


function p = setup_parameters(c, p_in)

  p = p_in;
  
  p.cmpt_and_wrt_sum_stats      = true;
  
  p.rescale_date                = c.Jan_2000;                 % p.rescale_date must be c.Jan_1953, c.Jan_2000, c.rescale_1987_last_mth_ave
  p.rescale_ave_mths            = 120;                        % relevant only if p.rescale_date = c.rescale_frst_5_yr_ave or c.rescale_frst_10_yr_ave
  
  if p.rescale_date == c.Jan_2000
    p.label_rescale_base          = "200001 Base";
  elseif p.rescale_date == c.Jan_1987
    p.label_rescale_base          = "198701 Base";
  elseif p.rescale_date == c.rescale_1987_last_mth_ave
    if p.last_mth == c.Sep_2019
      p.label_rescale_base        = "1987-2019 Ave Base";
    else
      p.label_rescale_base        = "1987-2020 Ave Base";
    end
  else
    disp('p.rescale_date must be Jan_2000, Jan_1987, or rescale_ave')
    pause
  end
  
  p.use_nominal_prices_and_rents     = false;
  
  p.PR_or_RP                         = c.price_to_rent;          % This is just for TS regression of dR on P/R and dP.
  p.ln_PR                            = false;

  p.wrt_just_US_average              = false;

  p.annualize_cc_change              = false;
  
  p.ave_chg_req_price_and_rent       = true;
  p.P_and_R_reqd_for_all_mthly_means = true;                      % If true, include only cities with both price and rent when computing all monthly means
  
  p.write_hetero_adj                 = false;
  p.wrt_raw_data                     = false;   % this is for diagnostics
																 
							   
  p.plot_city_real                   = true;   % *** true;
  p.city_plot_just_all               = true;
  p.plot_city_price_rent_ratio       = true;
							   
  p.plots_with_lines                 = true;
							   
  p.min_mths_partial_period          = 24;
							   
  p.min_mths_in_corrs_etc            = 24;
  p.num_autos                        = 12;
  p.num_autos_sum_stats              = 18;
  
end


function p = setup_output_file_etc(c, d, p)

  if p.use_nominal_prices_and_rents
    nominal_or_real = "nominal";
  else
    nominal_or_real = "real";
  end

  label_last = "Oct_2022";
  date = 20231030;

  if p.annualize_cc_change
    p.fid = fopen(sprintf('City_Sum_Stats_and_Plots_%s_annualized_%s_%8d.txt', nominal_or_real, label_last, date), 'wt');
  else
    p.fid = fopen(sprintf('City_Sum_Stats_and_Plots_%s_%s_%8d.txt', nominal_or_real, label_last, date), 'wt');
  end
    
  if p.SA_prices
    fprintf(p.fid, '\nPrices are seasonally adjusted.');
  else
    fprintf(p.fid, '\nPrices are not seasonally adjusted.');
  end
  
  if p.SA_rents
    fprintf(p.fid, '\nRents for All are seasonally adjusted; Seasonally adjusted data are not available for individual cities.');
  else
    fprintf(p.fid, '\nRents for All and individual cities are not seasonally adjusted.');
  end
  
  fprintf(p.fid, '\n\nRents are the Rent of Primary Residence component of the CPI.');
 
 
  if p.PR_or_RP == c.price_to_rent
    fprintf(p.fid, '\n\nThis run uses price/rent, not rent/price.');
  else
    fprintf(p.fid, '\n\nThis run uses rent/price, not price/rent.');
  end
  
  if p.ln_PR
    fprintf(p.fid, '\n\nThe ratios of price to rent (or rent to price) are in logs.');
  else
    fprintf(p.fid, '\n\nThe ratios of price to rent (or rent to price) are NOT in logs.');
  end    
  
  if p.rescale_date == c.rescale_1987_last_mth_ave
    if c.last_mth == c.Sep_2019
      fprintf(p.fid, '\n\nPrices, rents, price/rent, and rent/price are rescaled so averages for 1987-2019 are 1.0');
    elseif c.last_mth == c.Apr_2020
      fprintf(p.fid, '\n\nPrices, rents, price/rent, and rent/price are rescaled so averages for 1987-2020 are 1.0');
    elseif c.last_mth == c.Oct_2020
      fprintf(p.fid, '\n\nPrices, rents, price/rent, and rent/price are rescaled so averages for 198701-202010 are 1.0');
    else
      fprintf(p.fid, '\n\nPrices, rents, price/rent, and rent/price are rescaled so averages for 1987 to last mth are 1.0');
    end
  elseif p.rescale_date == c.rescale_frst_5_yr_ave
    fprintf(p.fid, '\n\nPrices and rents are rescaled so averages for 1987-1991 are 1.0');
  elseif p.rescale_date == c.rescale_frst_10_yr_ave
    fprintf(p.fid, '\n\nPrices and rents are rescaled so averages for 1987-1996 are 1.0');
  elseif 0 < p.rescale_date
    fprintf(p.fid, '\n\nPrices and rents are rescaled so values for %6d are 1.0', d.idate(p.rescale_date));
  end
 
end
  

function [price_filename, rent_filename] = get_input_file_names(c, p)

  if p.SA_prices
    price_filename = 'CS_Prices_SA_198701_202210.txt';
  else
    price_filename = 'CS NSA Prices 198701_202010.txt';
    disp('Non-seasonal CS Prices were not downloaded.')
    pause
  end
  
  rent_filename = 'CPI_Rent_of_Primary_Residence_NSA_198701_202211.txt';
end
  

function d = load_CPI_price_and_rent_data(price_filename, rent_filename, c, p)


  % Read CPI for January 1987 to July 2021
  %
  CPI = dlmread('CPI_Monthly_194701_202211.txt', '', [481, 0, 480 + c.last_mth, 2]);   % date,  NSA,  SA

  if p.CPI_SA     % true: real prices and rents  divided by CPI_SA
    d.CPI = CPI(:, 3);
  else
    d.CPI = CPI(:, 2);
  end

  d.idate   = CPI(:, 1);
  d.date    = floor(d.idate ./ 100) + (2 * mod(d.idate, 100) - 1) ./ 24;
  

  % setup price and rent data
  %
  d.price     = NaN(c.num_mths_ttl, c.num_cities_price);
  d.rent      = NaN(c.num_mths_ttl, c.num_cities_both);
  
  d.cityname  = reshape(blanks(c.num_cities_price * 13), c.num_cities_price, 13);
  d.cityabrev = reshape(blanks(c.num_cities_price * 5),  c.num_cities_price, 5);

  d.frst_ob   = NaN(c.num_cities_price, 3);
  d.last_ob   = NaN(c.num_cities_price, 3);

 
  % read price data  
  %
  f_in = fopen(price_filename, 'r');
  
  cityname_in  = textscan(f_in, [repmat('%s',1,c.num_cities_price), '%*[^\n]'], 1, 'HeaderLines', 34); 
  cityabrev_in = textscan(f_in, [repmat('%s',1,c.num_cities_price), '%*[^\n]'], 1); 
  frst_ob_in   = textscan(f_in, [repmat('%d',1,c.num_cities_price), '%*[^\n]'], 1);
  last_ob_in   = textscan(f_in, [repmat('%d',1,c.num_cities_price), '%*[^\n]'], 1);

  cityname  = reshape(blanks(c.num_cities_price * 13), c.num_cities_price, 13);
  cityabrev = reshape(blanks(c.num_cities_price * 5),  c.num_cities_price, 5);
  frst_ob   = cell2mat(frst_ob_in);
  last_ob   = cell2mat(last_ob_in);
  
  for row = 1:c.num_cities_price
    cityname (row, :) = strjust(sprintf('%13s', cell2mat(cityname_in {1, row})), 'left');
    cityabrev(row, :) = strjust(sprintf('%5s',  cell2mat(cityabrev_in{1, row})), 'left');
    frst_ob(row)      = loc_mth(frst_ob(row), c);
    last_ob(row)      = loc_mth(last_ob(row), c);
  end

  
  price_in = textscan(f_in, ['%d', repmat('%f',1,c.num_cities_price), '%*[^\n]'], c.num_mths_ttl, 'CollectOutput',1);
  price_idate = price_in{1};
  price       = price_in{2};

  loc_has_data = false(c.num_cities_price, 1);  
  
  for city = 1:c.num_cities_price
    loc = find_cityabrev_match(cityabrev(city, :), c);

    if loc == 0
      disp('Cant find city abreviation');
      disp(city);
      disp(cityabrev(city, :));
      pause;

    elseif c.city_loc_in_data(loc) ~= 0
      loc_has_data(c.city_loc_in_data(loc)) = true;
      
      d.price (:, c.city_loc_in_data(loc))          = price    (:, city);    % c.city_loc_in_data is where the city is in c.order of cities
      d.frst_ob  (c.city_loc_in_data(loc), c.price) = frst_ob  (city);       % frst_ob(:, 1) and last_ob(:, 1) is for price
      d.last_ob  (c.city_loc_in_data(loc), c.price) = min(last_ob(city), c.last_mth) ;
      d.cityname (c.city_loc_in_data(loc), :)       = cityname (city, :);
      d.cityabrev(c.city_loc_in_data(loc), :)       = cityabrev(city, :);
    end
  end
  fclose(f_in);
  
  % read rent data  
  %
  f_in = fopen(rent_filename, 'r');
  
  cityname_in  = textscan(f_in, [repmat('%s',1,c.num_cities_both + 1), '%*[^\n]'], 1, 'HeaderLines', 1); 
  cityabrev_in = textscan(f_in, [repmat('%s',1,c.num_cities_both + 1), '%*[^\n]'], 1); 
  frst_ob_in   = textscan(f_in, [repmat('%d',1,c.num_cities_both + 1), '%*[^\n]'], 1);
  last_ob_in   = textscan(f_in, [repmat('%d',1,c.num_cities_both + 1), '%*[^\n]'], 1);

  cityabrev = reshape(blanks(c.num_cities_both * 5), c.num_cities_both, 5);
  frst_ob   = cell2mat(frst_ob_in);
  last_ob   = cell2mat(last_ob_in);
  
  for row = 1:c.num_cities_both + 1
    cityname (row, :) = strjust(sprintf('%13s', cell2mat(cityname_in {1, row})), 'left');
    cityabrev(row, :) = strjust(sprintf('%5s',  cell2mat(cityabrev_in{1, row})), 'left');
    frst_ob  (row)    = loc_mth(frst_ob(row), c);
    last_ob  (row)    = min([loc_mth(last_ob(row), c), c.last_mth]);
  end
  
  rent_in    = textscan(f_in, ['%d', repmat('%f',1,c.num_cities_both + 1), '%*[^\n]'], c.num_mths_ttl, 'CollectOutput',1);
  rent_idate = rent_in{1};
  rent       = rent_in{2};
  

  Check_Alignment(d.idate, price_idate, rent_idate);
 
 
  % Load rents in columns defined by c.order_of_cities_stored
  %
  for city = 1:c.num_cities_both + 1
  
    if (p.SA_rents && strcmp(cityabrev(city, :),'SAAve')) || (~p.SA_rents && strcmp(cityabrev(city, :),'NSAAv'))
      loc = find_cityabrev_match('USAve', c);  % loc is the location of the city in c.cityname, which is also the value of the city's name, e.g. c.USAve = 1
      
    elseif ~strcmp(cityabrev(city, :),'SAAve') && ~strcmp(cityabrev(city, :),'NSAAv')    % want only one of SA Ave and NSA Ave; that one is selected in if statement above; skip the other
      loc = find_cityabrev_match(cityabrev(city, :), c);
    end

    if loc == 0
      disp('Cant find city abreviation');
      disp(city);
      disp(cityabrev(city, :));
      pause;

    elseif c.city_loc_in_data(loc) ~= 0 && c.city_loc_in_data(loc) <= c.num_cities_both
      d.rent(:, c.city_loc_in_data(loc))         = rent(:, city);
      d.frst_ob(c.city_loc_in_data(loc), c.rent) = frst_ob(city);    % frst_ob(:, 1) and last_ob(:, 1) is for price
      d.last_ob(c.city_loc_in_data(loc), c.rent) = last_ob(city);

      if ~loc_has_data(c.city_loc_in_data(loc))  % loc_has_data = t if price data available for city
        disp('City in rent file is in c.order_of_cities_stored but was not found in price file.')
        disp(cityabrev(city, :))
        pause;
        
        % Here because city has no price data
        %
        loc_has_data(c.city_loc_in_data(loc))    = true;
        d.cityname  (c.city_loc_in_data(loc), :) = cityname (city, :);
        d.cityabrev (c.city_loc_in_data(loc), :) = cityabrev(city, :);

      end
    end
  end
  fclose(f_in);


  % Set values of price or rent <= 0 to NaN
  %
  d.price(d.price <= 0) = NaN;
   d.rent(d.rent  <= 0) = NaN;


  % build header with city abbreviations
  %
  d.header_cities = reshape(blanks(c.num_cities_price * 8), c.num_cities_price, 8);
  for row = 1:c.num_cities_price
    d.header_cities(row, :) = sprintf('   %5s', strjust(d.cityabrev(row, :), 'right'));
  end      


  % Compute ratio of rent to price and price to rent, and the average values of the ratios (used in the const_x regression)
  %
  d.rent_price           = NaN(c.last_mth, c.num_cities_both);
  d.price_rent           = NaN(c.last_mth, c.num_cities_both);
  d.mean_rent_price_city =          NaN(1, c.num_cities_both);
  d.mean_price_rent_city =          NaN(1, c.num_cities_both);

  d.mean_rent_price_mth  = NaN(c.last_mth, 1);
  d.mean_price_rent_mth  = NaN(c.last_mth, 1);
  
  for city = 1:c.num_cities_both
    if 0 < min(d.frst_ob(city, c.price:c.rent))
      d.frst_ob(city, c.ratio) = max(d.frst_ob(city, c.price:c.rent));
      d.last_ob(city, c.ratio) = min(d.last_ob(city, c.price:c.rent));
      
      frst_mth = d.frst_ob(city, c.ratio);
      last_mth = d.last_ob(city, c.ratio);
      if frst_mth <= last_mth
        if p.ln_PR
          d.rent_price(frst_mth:last_mth, city) = log(d.rent (frst_mth:last_mth, city) ./ d.price(frst_mth:last_mth, city));
          d.price_rent(frst_mth:last_mth, city) = log(d.price(frst_mth:last_mth, city) ./ d.rent (frst_mth:last_mth, city));
        else
          d.rent_price(frst_mth:last_mth, city) = d.rent (frst_mth:last_mth, city) ./ d.price(frst_mth:last_mth, city);
          d.price_rent(frst_mth:last_mth, city) = d.price(frst_mth:last_mth, city) ./ d.rent (frst_mth:last_mth, city);
        end
        
        d.mean_price_rent_city(1, city) = mean(d.price_rent(frst_mth:last_mth, city));
        d.mean_rent_price_city(1, city) = mean(d.rent_price(frst_mth:last_mth, city));
      end
    end
  end    


  % standardize series...January 2000 is 1.0...note, this screws up -99.999, can only check whether value is negative
  %
  % In the rest of the project we typically rescale prices and rents so they are 1 on Jan 2000; p.rescale_date is usually Jan_2000
  % Here the base is usually the average value from 1987-2020. This choice affects only the FM regressions. 
  % 
  if 0 < p.rescale_date   % 0 < p.rescale_date means rescale so prices and rents are 1 on p.rescale_date
    for city = 1:c.num_cities_price
      if d.frst_ob(city, c.price) ~= 0 && d.frst_ob(city, c.price) <= p.rescale_date && p.rescale_date <= d.last_ob(city, c.price)
        d.price(d.frst_ob(city, c.price):d.last_ob(city, c.price), city) = d.price(d.frst_ob(city, c.price):d.last_ob(city, c.price), city) ./ d.price(p.rescale_date, city);
      else
        d.frst_ob(city, c.price) = 0;
      end
    end

    for city = 1:c.num_cities_both
      if d.frst_ob(city, c.rent) ~= 0 && d.frst_ob(city, c.rent) <= p.rescale_date && p.rescale_date <= d.last_ob(city, c.rent)
        d.rent(d.frst_ob(city, c.rent):d.last_ob(city, c.rent), city) = d.rent(d.frst_ob(city, c.rent):d.last_ob(city, c.rent), city) ./ d.rent(p.rescale_date, city);
      else
        d.frst_ob(city, c.rent) = 0;
      end
    end    
    
    for city = 1:c.num_cities_both
      if d.frst_ob(city, c.ratio) ~= 0 && d.frst_ob(city, c.ratio) <= p.rescale_date && p.rescale_date <= d.last_ob(city, c.ratio)
        d.rent_price(d.frst_ob(city, c.ratio):d.last_ob(city, c.ratio), city) = d.rent_price(d.frst_ob(city, c.ratio):d.last_ob(city, c.ratio), city) ./ d.rent_price(p.rescale_date, city);
        d.price_rent(d.frst_ob(city, c.ratio):d.last_ob(city, c.ratio), city) = d.price_rent(d.frst_ob(city, c.ratio):d.last_ob(city, c.ratio), city) ./ d.price_rent(p.rescale_date, city);
      else
        d.frst_ob(city, c.ratio) = 0;
      end
    end    
    
    d.CPI = d.CPI ./ d.CPI(p.rescale_date, :);

  else
  
    for city = 1:c.num_cities_price
      if d.frst_ob(city, c.price) ~= 0
        city_mean_price = mean(d.price(d.frst_ob(city, c.price):d.last_ob(city, c.price), city));
        d.price(d.frst_ob(city, c.price):d.last_ob(city, c.price), city) = d.price(d.frst_ob(city, c.price):d.last_ob(city, c.price), city) ./ city_mean_price;
      end
      
      if d.frst_ob(city, c.rent) ~= 0
        city_mean_rent = mean(d.rent(d.frst_ob(city, c.rent):d.last_ob(city, c.rent), city));
        d.rent(d.frst_ob(city, c.rent):d.last_ob(city, c.rent), city) = d.rent(d.frst_ob(city, c.rent):d.last_ob(city, c.rent), city) ./ city_mean_rent;
      end

      if d.frst_ob(city, c.ratio) ~= 0
        frst_mth = d.frst_ob(city, c.ratio);
        last_mth = d.last_ob(city, c.ratio);
        if frst_mth <= last_mth
          d.price_rent(frst_mth:last_mth, city) = d.price_rent(frst_mth:last_mth, city) ./ d.mean_price_rent_city(1, city);
          d.rent_price(frst_mth:last_mth, city) = d.rent_price(frst_mth:last_mth, city) ./ d.mean_rent_price_city(1, city);
        end
      end
    end

    d.CPI = d.CPI ./ mean(d.CPI(c.Jan_1987:end));
  end
  

  % compute real prices and rents
  %
  d.real_price = d.price ./ d.CPI;
  d.real_rent  = d.rent  ./ d.CPI;

  if p.PR_or_RP == c.price_to_rent
    d.PR = d.price_rent;
  else
    d.PR = d.rent_price;
  end
end


function loc = find_cityabrev_match(cityabrev, c)
  
  for loc = 1:size(c.cityabrev, 1)
    if strcmp(c.cityabrev(loc, :), cityabrev)
      return
    end
  end
  loc = 0;
end


%-----------------------------------------------------------------------------------------
%                                                                                        *
%                                   Functions to Plot Data                               *
%                                                                                        *
% ****************************************************************************************
%
function cntrl_plot_city_real(city, c, d, p)


  label = " ";
  price      = d.real_price(:, city);
  rent       = d.real_rent (:, city);
  price_rent = d.price_rent(:, city);
  rent_price = d.rent_price(:, city);
  
  frst_obs = [d.frst_ob(city, c.price), d.frst_ob(city, c.rent), d.frst_ob(city, c.ratio)];
  last_obs = [d.last_ob(city, c.price), d.last_ob(city, c.rent), d.last_ob(city, c.ratio)];
  
  
  % plot prices, rents, and ratios
  %
  filename = sprintf('Fig 1 Prices and Rents, %s, %s', p.label_rescale_base, strtrim(c.full_city_name(city, 1)));
  
  plot_three_real_variables(d.date, frst_obs, last_obs, price, rent, price_rent, filename, label, c, p);
  
end


function plot_three_real_variables(date, frst_obs, last_obs, x1, x2, x3, filename, label, c, p)

  % frst_obs and last_obs are for prices and rents
  % x1 is prices, x2 is rents, x3 is price/rent
  %
  hold off

  plot(date(frst_obs(c.price):last_obs(c.price)), x1(frst_obs(c.price):last_obs(c.price)), '-k', "LineWidth", 1.0);
  hold on
  plot(date(frst_obs(c.rent):last_obs(c.rent)),   x2(frst_obs(c.rent):last_obs(c.rent)),  '--b', "LineWidth", 1.25);
  plot(date(max(frst_obs):min(last_obs)),         x3(max(frst_obs):min(last_obs)),   ':r', "LineWidth", 2.0);
  
  
    x0 = 0.1;
    y0 = 0.1;
    width  = 7.0;  %6.5
    height = 2.75;  %2.6
    set(gcf,'unit', 'inches', 'position', [x0, y0, width, height])
    
    xlabel('Date', 'FontSize', 30, 'FontWeight', 'bold', 'Color', 'black')
    ylabel('Index', 'FontSize', 30, 'FontWeight', 'bold', 'Color', 'black')


    if p.city_plot_just_all
      min_y = 0.8;
      max_y = 1.8;
      yticks([0.8:0.2:1.8])
    else
      min_y = 0.5;
      max_y = 2.5;
    end
    
    xlim([min(date) max(date)])
    ylim([min_y max_y]);
    
    if p.plots_with_lines
      coord_x = [date(c.real_key_obs([2, 4, 6]))', min(date), min(date); date(c.real_key_obs([2, 4, 6]))', max(date), max(date)];
      coord_y = [min_y * ones(1, 3), 0, 1; max_y * ones(1, 3), 0, 1];
    else
      coord_x = [min(date), max(date)];
      coord_y = [0, 0];
    end
    line(coord_x, coord_y, 'Color', 'black')

    if ~p.city_plot_just_all
      text(min(date) + 1, 2.35, label, 'FontSize', 30, 'FontWeight', 'bold', 'Color', 'black');
%    text(min(date) + 1, 1.35, label, 'FontSize', 30, 'FontWeight', 'bold', 'Color', 'black');
    end
      
    ax = gca;
    ax.FontSize = 24;
    ax.FontWeight = 'bold';

    savefig(filename);
  
    hold off
end
      
%-----------------------------------------------------------------------------------------
%                                                                                        *
%                       Utility Functions used in Old and/or New                         *
%                                                                                        *
% ****************************************************************************************
%
% OK
%
function loc = loc_mth(date, c)   % date is yyyymm, like 201012   loc is the location relative to the c.base_date, like 198701
  
  cal_yr   = floor(date ./ 100);
  cal_mth  = mod(date, 100);
  base_yr  = floor(c.base_date ./ 100);
  base_mth = mod(c.base_date, 100);

  loc      = 12 * (cal_yr - base_yr) + (cal_mth - base_mth) + 1;
end
  

% OK
%
function write_raw_data(c, d, p)

  fprintf(p.fid, '\n\n\nCS Prices');
  fprintf(p.fid, '\n\n      '); 
  for city = 1:c.num_cities_price
    fprintf(p.fid, ' %s', c.name_7(c.order_price_stored(city), :)); 
  end

  for j = 1:size(d.price, 1)
    fprintf(p.fid, '\n%6d', d.idate(j, 1));
    fprintf(p.fid, '%8.3f', d.price(j, :));
  end
  
  fprintf(p.fid, '\n\n\nCPI Rents and CPI');
  fprintf(p.fid, '\n\n      '); 
  
  for city = 1:c.num_cities_both
    fprintf(p.fid, ' %s', c.name_7(c.order_rent_stored(city), :)); 
  end
  fprintf(p.fid, '  CPI_SA CPI_NSA');
  
  for j = 1:size(d.rent, 1)
    fprintf(p.fid, '\n%6d', d.idate(j, 1));
    fprintf(p.fid, '%8.3f', d.rent(j, :)); fprintf(p.fid, '%8.3f', d.CPI_SA(j), d.CPI_NSA(j));
  end

end

% OK
%
function Check_Alignment(idate, price_idate, rent_idate)

  % Make sure files are aligned
  %
  if size(rent_idate, 1) ~= size(price_idate, 1)
    disp('sizes of rent_idate and price_idate do not match')
    size(rent_idate), size(price_idate)
    pause
  end
  
  if size(rent_idate, 1) ~= size(idate, 1)
    disp('sizes of rent_idate and cpi idate do not match')
    size(rent_idate), size(idate)
    pause
  end
  
  for j = 1:size(idate, 1)
    if rent_idate(j, 1) ~= price_idate(j, 1)
      disp('rent_idate and price_idate misaligned')
      rent_idate(j, 1), price_idate(j, 1)
      pause
    elseif idate(j, 1) ~= rent_idate(j, 1)
      disp('rent_idate and d.idate misaligned')
      rent_idate(j, 1), idate(j, 1)
      pause
    end
  end
end  


% OK_isnan
%
function [frst_obs, last_obs] = get_frst_last_obs(idate, x, c)

  IDX = uint32(1:size(x, 1));

  frst_obs = NaN(size(x, 2), 1);
  last_obs = NaN(size(x, 2), 1);
  
  for col = 1:size(x, 2)
    if any(~isnan(x(:, col)))
      frst_obs(col) = loc_mth(idate(min(IDX(~isnan(x(:, col))))), c);
      last_obs(col) = min(loc_mth(idate(max(IDX(~isnan(x(:, col))))), c), c.last_mth);
    end
  end
end




